*SANTOS, ALMEIDA & MOREIRA (BPSR, 2019)

*ESTIMATION PROCEDURES

*Define dependent variable (y)
local y= "gtot"

*Define vectors of independent variables (x) and number of vectors (I)
local x1= "idepre vetdis finopec ncommc trdope curacc margin unemp gdpkpl"
local x2= "idepre vetdis finopec ncommc trdope curacc votpre unemp gdpkpl"
local x3= "idepre vetdis finopec trdope curacc votpre unemp gdpkpl vet_fope"
local x4= "idepre vetdis ncommc trdope curacc votpre unemp gdpkpl vet_ncom"
local x= "idepre vetdis finope ncomm trdope curacc votpre margin unemp gdpkpl"
local I=4

*Define number of imputed datasets
local M= 5

quietly {

est drop _all
mat drop _all

*Create matrix of test results
mat TEST= (0,0,0,0,0,0,0,0)
mat coln TEST= MEsp DSet THet IHet IRE Hetk CCor SCor

local i=1				//model specification index
while `i'<=`I' {

local m=1				//data file index
while `m'<=`M' {

*Open data file (user-defined)
use C:\Users\r1702830\Documents\midtset`m'.dta, clear

*Declare data as panel
xtset id year


*PART 1: PRELIMINARY TESTS

*F test of temporal heterogeneity (thet)
xtreg `y' `x`i'' i.year
testparm i.year
scalar thet= r(p)

*Breusch-Pagan test of cross-sectional heterogeneity (ihet)
xtreg `y' `x`i'', re
est store re
xttest0
scalar ihet= r(p)

*Mundlak test of random effects (ire)
if ihet< .20 {
mundlak `y' `x`i'', keep
testparm mean_*
scalar ire= r(p)
}
else {
scalar ire= .
}
*

*Likelihood-ratio test of panel heteroskedasticity
xtgls `y' `x`i'', igls panels(heteroskedastic)
est store hetk
xtgls `y' `x`i'', igls
est store iid
local df= e(N_g)-1
lrtest hetk iid, df(`df')
scalar hetk= r(p)

*Pesaran's test of contemporaneous correlation
est restore iid
predict yhat if e(sample)
gen res= yhat-`y'
xtcd res, resid
scalar ccor= el(r(pesaran),1,1)
drop yhat res

*Wooldridge's test of panel serial correlation
xtserial `y' `x`i''
scalar scor= r(p)

*Fill in matrix with test results 
mat TEST`i'`m'= (`i',`m',thet,ihet,ire,hetk,ccor,scor)
mat TEST= TEST\TEST`i'`m'


*Levin-Lin-Chu test of panel unit-root
if `i'==1 {
foreach var of varlist gtot gpro ghea gedu {
xtunitroot llc `var'
gen series_`var'= r(p_tds)
}
*

*Create matrix with test results
tabstat series_*, s(mean) va(16) save
mat UROOT`m'= r(StatTotal)'
}
*

*Variance inflation factors for pooled (VIF) and least squares dummy variable (VIFD) models
if `m'==1 {
reg `y' `x`i''
local nvar= e(df_m)
estat vif

*Create matrix with results
mat VIF`i'= (0)
mat rown VIF`i'= fstrow
forvalues v= 1/`nvar' {
mat M`v'= r(vif_`v')
mat rown M`v'= `r(name_`v')'
mat VIF`i'= VIF`i'\M`v'
}
*

reg `y' `x`i'' i.id
local nvar= e(df_m)
estat vif

*Create matrix with results
mat VIFD`i'= (0)
mat rown VIFD`i'= fstrow
forvalues v= 1/`nvar' {
mat M`v'= r(vif_`v')
mat rown M`v'= `r(name_`v')'
mat VIFD`i'= VIFD`i'\M`v'
}
*

*Adjust matrices
mat VIF`i'= VIF`i'[2...,1...]
mat VIFD`i'= VIFD`i'[2...,1...]
mat coln VIF`i'= DS`m'
mat coln VIFD`i'= DS`m'
}
*

*PART 2: MODEL ESTIMATION

*Estimate static random effects w/ Driscoll-Kraay std errors (DKRE)
xtscc `y' `x`i'', lag(1) re
est store DKRE`i'`m'

*Create matrices with estimates (BRE) and variance (VRE)
mat BRE= e(b)'
mat VRE= vecdiag(e(V))
mat VRE= VRE'

*Calculate adjusted R2
scalar nre= e(N)
scalar rsqre= 1-((1-e(r2_o))*(e(N)-1)/(e(N)-e(df_m)-1))


*Estimate dynamic random effects w/ Driscoll-Kraay std errors (LDV)
xtscc `y' `x`i'' `y'lag, re
est store LDV`i'`m'

*Create matrices with estimates (BLDV) and variance (VLDV)
mat A= e(b)'
scalar nrow= rowsof(A)
scalar bx= nrow-2
mat B= A[1..bx,1]
mat C= A[nrow,1]
mat rown C= _cons
mat BLDV= B\C
scalar lag= nrow-1
scalar blag= el(A,lag,1)
mat EV= vecdiag(e(V))
mat A= EV'
mat B= A[1..bx,1]
mat C= A[nrow,1]
mat rown C= _cons
mat VLDV= B\C
scalar vlag= el(A,lag,1)
	
*Calculate adjusted R2
scalar nldv= e(N)
scalar rsqldv= 1-((1-e(r2_o))*(e(N)-1)/(e(N)-e(df_m)-1))


*Estimate fixed effects w/ Driscoll-Kraay std errors
xtscc `y' `x`i'', lag(1) fe
	
*Create matrices with estimates (BFE) and variance (VFE)
est store DKFE`i'`m'
mat BFE= e(b)'
mat VFE= vecdiag(e(V))
mat VFE= VFE'
		
*Calculate adjusted R2
scalar nfe= e(N)
scalar rsqfe= 1-((1-e(r2_w))*(e(N)-1)/(e(N)-e(df_m)-1))


*Estimate fixed effects with panel corrected std errors
xtpcse `y' `x`i'' i.id, c(ar1)
est store PCSE`i'`m'
		
*Create matrices with estimates (BPCSE) and variance (VPCSE)
mat A= e(b)'
scalar nrow= rowsof(A)
tab id if e(sample)
scalar nid= r(r)
scalar bx= nrow-nid-1
mat B= A[1..bx,1]
mat C= A[nrow,1]
mat rown C= _cons
mat BPCSE= B\C
mat EV= vecdiag(e(V))
mat A= EV'
mat B= A[1..bx,1]
mat C= A[nrow,1]
mat rown C= _cons
mat VPCSE= B\C
		
*Calculate adjusted R2
scalar npcse= e(N)
scalar rsqpcse= 1-((1-e(r2))*(e(N)-1)/(e(N)-e(df_m)-1))


*Create matrices with results from all models
mat B= BRE,BLDV,BFE,BPCSE
mat V= VRE,VLDV,VFE,VPCSE
mat BLAG= (0,blag,0,0)
mat VLAG= (0,vlag,0,0)
mat rown BLAG= `y'lag1
mat rown VLAG= `y'lag1
mat B`i'`m'= B\BLAG
mat V`i'`m'= V\VLAG
mat N`i'`m'= (nre,nldv,nfe,npcse)
mat RSQ`i'`m'= (rsqre,rsqldv,rsqfe,rsqpcse)

local m= `m'+1
}
*

*PART 3: POOLED MULTIPLE-IMPUTATION ESTIMATES AND MEASURES (Schafer, 1997, formulae 4.21-4.26)

*Create matrix with MI estimates (Q)
mat Q`i'= (B`i'1+B`i'2+B`i'3+B`i'4+B`i'5)/5

*Create matrix with MI variances (U)
mat WIU`i'= (V`i'1+V`i'2+V`i'3+V`i'4+V`i'5)/5				//within-imputation variance

forvalues m= 1/5 {
mat D`i'`m'= B`i'`m'-Q`i'
mat DSQ`i'`m'= hadamard(D`i'`m',D`i'`m')
}
mat BIU`i'= (DSQ`i'1+DSQ`i'2+DSQ`i'3+DSQ`i'4+DSQ`i'5)/4		//between-imputation variance
mat U`i'= WIU`i'+(1+1/5)*BIU`i'								//total variance

*Create matrices with MI std errors (SE), t statistics (T), degrees of freedom (DF), and p-values (PV)
mat SE`i'= U`i'
mat T`i'= U`i'
mat DF`i'= U`i'
mat PV`i'= U`i'
forvalues v= 1/`= rowsof(U`i')' {
forvalues vv= 1/`= colsof(U`i')' {
mat SE`i'[`v',`vv']= sqrt(U`i'[`v',`vv'])
mat T`i'[`v',`vv']= abs(Q`i'[`v',`vv']/SE`i'[`v',`vv'])
mat DF`i'[`v',`vv']= round(4*(1+WIU`i'[`v',`vv']/((1+1/5)*BIU`i'[`v',`vv']))^2)
mat PV`i'[`v',`vv']= ttail(DF`i'[`v',`vv'],T`i'[`v',`vv'])*2
}
} 
*
*Create matrices with sample sizes (N) and adjusted R2 (RSQ)
mat N`i'= (N`i'1+N`i'2+N`i'3+N`i'4+N`i'5)/5
mat RSQ`i'= (RSQ`i'1+RSQ`i'2+RSQ`i'3+RSQ`i'4+RSQ`i'5)/5
mat rown N`i'= N
mat rown RSQ`i'= R2ajt

*Create matrices with all estimates (EST)
mat ESTB`i'= Q`i'\N`i'\RSQ`i'
mat ESTSE`i'= SE`i'\N`i'\RSQ`i'
mat ESTPV`i'= PV`i'\N`i'\RSQ`i'
mat EST`i'= ESTB`i',ESTSE`i',ESTPV`i'

if `= colsof(EST`i')'== 6 {
mat coln EST`i'= DKRE_b LDV_b DKRE_se LDV_se DKRE_pv LDV_pv
}
else {
mat coln EST`i'= DKRE_b LDV_b DKFE_b PCSE_b DKRE_se LDV_se DKFE_se PCSE_se DKRE_pv LDV_pv DKFE_pv PCSE_pv
}
*

local i= `i'+1
}
*

*Create matrix with all unit-root test results
mat UROOT= UROOT1,UROOT2,UROOT3,UROOT4,UROOT5
mat coln UROOT= DS1 DS2 DS3 DS4 DS5

*Adjust matrix
mat TEST= TEST[2...,1...]
}
*

forvalues i= 1/`I' {
if `i'==1 {
matlist TEST, title("Test results") format(%8.0g) names(columns)
di ""
di "MEsp= Model specification"
di "DSet= Multiple-imputation dataset"
di "THet= Temporal heterogeneity F test"
di "IHet= Breusch-Pagan test of cross-sectional heterogeneity"
di "IRE= Mundlak test of random (vs. fixed) effects"
di "Hetk= Panel heteroskedasticity LR test"
di "CCor= Pesaran's test of contemporaneous correlation (Ho: 0<P<1)"
di "SCor= Wooldridges test of panel serial correlation"
di ""
di ""
matlist UROOT, title("Levin-Lin-Chu test of panel unit-root") format(%8.0g)
di ""
di "DS= Dataset"
di ""
di ""
}
*
di "Model specification: "`i'
matlist EST`i',  format(%8.0g)
di ""
di "DKRE= random effects w/ Driscoll-Kraay std errors and AR1"
di "LDV= lagged dep var random effects w/ Driscoll-Kraay std errors"
di "DKFE= fixed effects w/ Driscoll-Kraay std errors and AR1"
di "PCSE= panel corrected standard errors w/ AR1"
di "b= beta; se= std error; pv= p-value"
di ""
di ""
}
*

forvalues i= 1/`I' {
di "Model specification: "`i'
matlist VIF`i', title("VIFs: without country dummies") format(%8.0g)
di ""
*matlist VIFD`i', title("VIFs: with country dummies") format(%8.0g)
*di ""
}
*

